CONTENTS

1. Introduction

Recent advances in next-generation sequencing, microarrays and mass spectrometry for omics data production have enabled the generation and collection of different modalities of high-dimensional molecular data1. Clustering multi-omic data has the potential to reveal further systems-level insights, but raises computational and biological challenges2. This vignette aims to show how to use the package named MOVICS to perform multi-omics integrative clustering and visualization for cancer subtyping researches. This R package provides a unified interface for 10 state-of-the-art multi-omics clustering algorithms, and standardizes the output for each algorithm so as to form a pipeline for downstream analyses. Ten algorithms are CIMLR, iClusterBayes, MoCluster, COCA, ConsensusClustering, IntNMF, LRAcluster, NEMO, PINSPlus, and SNF where the former three methods can also perform the process of feature selection. For cancer subtyping studies, MOVICS also integrates many downstream analyses for further subtype characterization and creates editable publication-quality illustrations (see more details below). Please note that MOVICS currently supports up to 6 omics data for jointly clustering and users must provide at least 2 omics datastes as input. Ok then, let us make our hands dirty to figure out how MOVICS works.

2. Installation

It is essential that you have R 3.6 or above already installed on your computer or server. MOVICS is a pipeline that utilises many other R packages that are currently available from CRAN and Bioconductor. For all of the steps of the pipeline to work, make sure that you have upgraded Bioconductor to newest version (BiocManager v3.10). After you have R and Bioconductor installed properly, you can start to install MOVICS. The easiest way to install it is by typing the following code into your R session:

When you are installing MOVICS, you may encounter some errors saying that some packages are not installed. These errors are caused by recursively depending on R packages, so if one package was not installed properly on your computer, MOVICS would fail. To solve these errors, you simply need to check those error messages, find out which packages required are missing, then install it with command BiocManager::install("YourErrorPackage") or install.packages("YourErrorPackage") directly. After that, retry installing MOVICS, it may take several times, but eventually it should work. After installation, you should be able to load the MOVICS package in your R session:

3. Real Data Scenario

This package contains two preprocessed real datasets of breast cancer. One dataset is brca.tcga.rda which is a list that includes 643 samples with four complete omics data types of breast cancer retrieved from TCGA-BRCA cohort3 (i.e., mRNA expression, lncRNA expression, DNA methylation profiling and somatic mutation matrix), and corresponding clinicopathological information (i.e., age, pathological stage, PAM50 subtype, vital status and overall suvival time); such data list also contains corresponding RNA-Seq raw count table and Fragments Per Kilobase Million (FPKM) data in order to test the functions for downstream analyses (e.g., differential expression analysis, drug sensitivity analysis, etc). The other one, brca.yau.rda, is an external validation dataset which contains gene expression profilings and clinicopathological information that downloaded from BRCA-YAU cohort4 with 682 samples (one sample without annotation of PAM50 subtype was removed), which can be used to test the predictive functions available in MOVICS. These two datasets can be loaded like below:

Since in most cases, multi-omics clustering procedure is rather time-consuming, the TCGA-BRCA dataset was threfore first pre-processed to extract top 500 mRNAs, 500 lncRNA, 1,000 promoter CGI probes/genes with high variation using statistics of median absolute deviation (MAD), and 30 genes that mutated in at least 3% of the entire cohort. All these features are used in the following clustering analyses.

4. MOVICS Pipeline

4.1 Pipeline Introduction

MOVICS Pipeline diagram above outlines the concept for this package and such pipeline comprises separate functions that can be run individually (see datails about argument list by typing ?function_name() in R session). In above graphical abstract, all functions of MOVICS are included, which are categorized into three modules using three colors:

  • GET Module: get subtypes through multi-omics integrative clustering
    • getElites(): get elites which are those features that pass the filtering procedure and are used for analyses
    • getClustNum(): get optimal cluster number by calculating clustering prediction index (CPI) and Gap-statistics
    • get%algorithm_name%(): get results from one specific multi-omics integrative clustering algorithm with detailed parameters
    • getMOIC(): get a list of results from multiple multi-omics integrative clustering algorithm with parameters by default
    • getConsensusMOIC(): get a consensus matrix that indicates the clustering robustness across different clustering algorithms
    • getStdiz(): get a standardized data for generating comprehensive multi-omics heatmap
    • getMoHeatmap(): get a comprehensive multi-omics heatmap based on clustering results

  • COMP Module: compare subtypes from multiple perspectives
    • compSurv(): compare survival outcome and generate a Kalan-Meier curve with pair-wise comparsion if possible
    • compClinvar(): compare and summarize clinical features among different identified subtypes
    • compMut(): compare mutational frequency and generate an oncoprint with significant mutations
    • compTMB(): compare total mutation burden among subtypes and generate distribution of Transitions and Transversions
    • compFGA(): compare fraction genome altered among subtypes and generate a barplot for distribution comparison
    • compDrugsen(): compare estimated half maximal inhibitory concentration (\(IC_{50}\)) for drug sensitity and generate a boxviolin for distribution comparison
    • compAgree(): compare agreement of current subtypes with other pre-existed classifications and generate an alluvial diagram and an agreement barplot

  • RUN Module: run marker identification and verify subtypes
    • runDEA(): run differential expression analysis with three popular methods for choosing, including edgeR, DESeq2, and limma
    • runMarker(): run biomarker identification to determine uniquely and significanly differential expressed genes for each subtype
    • runGSEA(): run gene set enrichment analysis (GSEA), calculate activity of functional pathways and generate a pathway-specific heatmap
    • runNTP(): run nearest template prediction based on identified biomarkers to evaluate subtypes in external cohorts

4.2 Steps of Pipeline

Basically, the above three connected modules explain the workflow of this R package. MOVICS first identifies the cancer subtype (CS) by using one or multiple clustering algorithms; if multiple clustering algorithms are specified, it is highly recommended to perform a consensus clustering based on different subtyping results in order to derive stable and robust subtypes. Second, after having subtypes it is natural to exploit the heterogeneity of subtypes from as many angles as possible. Third, each subtype should have a list of subtype-specific biomarkers for reproducing such subtype in external cohorts. Therefore, let us follow this analytic pipeline and step into each module and each function in MOVICS.

4.2.1 GET Module

1) get data from example files

Data used for this vignette should be first extracted from the list of brca.tcga, including 4 types of multi-omics data. Except for omics data, this list also contains RNA-Seq raw count, FPKM matrix, MAF data, segmented copy number, clinical and survival information for downstream analyses. Notably, all these omics data used for clustering share exactly the same samples with exactly the same order, and please make sure of this when using MOVICS on your own data.

2) get elites by reducing data dimension

Although all these omics data have been already processed (filtered from the original dataset), I am still pleased to show you how to use getElites() function to filter out features that meet some stringent requirements, and those features that are preserved in this procedure are considered elites by MOVICS. Four filtering methods are provided here, namely mad for median absolute deviation, sd for standard deviation, cox for univariate Cox proportional hazards regression, and none for nothing but is really necessary for binary omics data. This function also handles missing values coded in NA by removing them directly or imputing them by \(k\) nearest neighbors using a Euclidean metric through argument of na.action. Let me show you how to use getElites() below:

# scenario 1: considering we are dealing with an expression data that have 2 rows with NA values
tmp       <- brca.tcga$mRNA.expr # get expression data
dim(tmp) # check data dimension
#> [1] 500 643
tmp[1,1]  <- tmp[2,2] <- NA # set 2 rows with NA values
tmp[1:3,1:3] # check data
#>         BRCA-A03L-01A BRCA-A04R-01A BRCA-A075-01A
#> SCGB2A2            NA          1.42          7.24
#> SCGB1D2         10.11            NA          5.88
#> PIP              4.54          2.59          4.35
elite.tmp <- getElites(dat       = tmp,
                       method    = "mad",
                       na.action = "rm", # NA values will be removed
                       elite.pct = 1) # elite.pct equals to 1 means all (100%) features after NA removal will be selected even using mad method
#> --2 features with NA values are removed.
#> missing elite.num then use elite.pct
dim(elite.tmp$elite.dat) # check dimension again and see that we have removed 2 rows with NA data
#> [1] 498 643

elite.tmp <- getElites(dat       = tmp,
                       method    = "mad",
                       na.action = "impute", # NA values will be imputed
                       elite.pct = 1) 
#> missing elite.num then use elite.pct
dim(elite.tmp$elite.dat) # all data kept
#> [1] 500 643
elite.tmp$elite.dat[1:3,1:3] # NA values have been imputed 
#>         BRCA-A03L-01A BRCA-A04R-01A BRCA-A075-01A
#> SCGB2A2         6.867         1.420          7.24
#> SCGB1D2        10.110         4.739          5.88
#> PIP             4.540         2.590          4.35

# scenario 2: considering we are dealing with continous data and use mad or sd to select elites
tmp       <- brca.tcga$mRNA.expr # get expression data with 500 features
elite.tmp <- getElites(dat       = tmp,
                       method    = "mad",
                       elite.pct = 0.1) # this time only top 10% features with high mad values are kept
#> missing elite.num then use elite.pct
dim(elite.tmp$elite.dat) # get 50 elite left
#> [1]  50 643

elite.tmp <- getElites(dat       = tmp,
                       method    = "sd",
                       elite.num = 100, # this time only top 100 features with high sd values are kept
                       elite.pct = 0.1) # this time elite.pct arugument will be disabled because elite.num has been already indicated.
#> elite.num has been provided then discards elite.pct.
dim(elite.tmp$elite.dat) # get 100 elites left
#> [1] 100 643

# scenario 3: considering we are dealing with data and use cox to select elite
tmp       <- brca.tcga$mRNA.expr # get expression data 
elite.tmp <- getElites(dat       = tmp,
                       method    = "cox",
                       surv.info = surv.info, # must provide surivival information with 'futime' and 'fustat'
                       p.cutoff  = 0.05,
                       elite.num = 100) # this time elite.num arugument will be disabled because cox method refers to p.cutoff to select elites
#> --all sample matched between omics matrix and survival data.
#> 5% 10% 15% 20% 25% 30% 35% 40% 45% 50% 55% 60% 65% 70% 75% 80% 85% 90% 95% 100%
dim(elite.tmp$elite.dat) # get 125 elites
#> [1] 125 643
table(elite.tmp$unicox$pvalue < 0.05) # 125 genes have nominal pvalue < 0.05 in univariate Cox regression
#> 
#> FALSE  TRUE 
#>   375   125

tmp       <- brca.tcga$mut.status # get mutation data 
elite.tmp <- getElites(dat       = tmp,
                       method    = "cox",
                       surv.info = surv.info, # must provide surivival information with 'futime' and 'fustat'
                       p.cutoff  = 0.05,
                       elite.num = 100) # this time elite.num arugument will be disabled because cox method refers to p.cutoff to select elites
#> --all sample matched between omics matrix and survival data.
#> 7% 13% 20% 27% 33% 40% 47% 53% 60% 67% 73% 80% 87% 93% 100%
dim(elite.tmp$elite.dat) # get 3 elites
#> [1]   3 643
table(elite.tmp$unicox$pvalue < 0.05) # 3 mutations have nominal pvalue < 0.05
#> 
#> FALSE  TRUE 
#>    27     3

# scenario 4: considering we are dealing with mutation data without using cox to select elites
tmp       <- brca.tcga$mut.status # get mutation data 
rowSums(tmp) # check mutation frequency
#> PIK3CA   TP53    TTN   CDH1  GATA3   MLL3  MUC16 MAP3K1  SYNE1  MUC12 
#>    208    186    111     83     58     49     48     38     33     32 
#>    DMD  NCOR1    FLG   PTEN   RYR2  USH2A  SPTA1 MAP2K4  MUC5B    NEB 
#>     31     31     30     29     27     27     25     25     24     24 
#>   SPEN  MACF1   RYR3    DST  HUWE1  HMCN1  CSMD1  OBSCN   APOB  SYNE2 
#>     23     23     23     22     22     22     21     21     21     21
elite.tmp <- getElites(dat       = tmp,
                       method    = "none", # must set as 'none'
                       elite.num = 80, # note: in this scenario elite.num refer to frequency of mutation
                       elite.pct = 0.1) # discard because elite.num has been already indicated
#> --method of none only supports binary omics data (e.g., somatic mutation matrix), and in this manner, elite.pct and elite.num are used to cut frequency.
#> elite.num has been provided then discards elite.pct.
rowSums(elite.tmp$elite.dat) # only genes that are mutated in over than 80 samples are kept as elites
#> PIK3CA   TP53    TTN   CDH1 
#>    208    186    111     83

elite.tmp <- getElites(dat       = tmp,
                       method    = "none", # must set as 'none'
                       elite.pct = 0.2) # note: in this scenario elite.pct refer to frequency of mutation / sample size
#> --method of none only supports binary omics data (e.g., somatic mutation matrix), and in this manner, elite.pct and elite.num are used to cut frequency.
#> missing elite.num then use elite.pct
rowSums(elite.tmp$elite.dat) # only genes that are mutated in over than 0.2*643=128.6 samples are kept as elites
#> PIK3CA   TP53 
#>    208    186

# get mo.data list just like below (not run)
# mo.data <- list(omics1 = elite.tmp$elite.dat,
#                 omics2 = ...)

Now I think I have made this clear for you concerning how to reduce data dimension for clustering analysis. Since the mo.data has been already prepared, I am going to stick on this data list, and take you further to MOVICS.

3) get optimal number for clustering

The most important parameter to estimate in any clustering study is the optimum number of clusters \(k\) for the data, where \(k\) needs to be small enough to reduce noise but large enough to retain important information. Herein MOVICS refers to CPI5 and Gaps-statistics6 to estimate the number of clusters by using getClustNum() function.

Figure 1. Identification of optimal cluster number by calculating Cluster Prediction Index (blue line) and Gaps-statistics (red line) in TCGA-BRCA cohort.

Figure 1. Identification of optimal cluster number by calculating Cluster Prediction Index (blue line) and Gaps-statistics (red line) in TCGA-BRCA cohort.

The above estimation of clustering number gives an arbitrary \(k\) of 2. However, the popular PAM50 classifier for breast cancer has 5 classifications and if taking a close look at the descriptive figure, it can be noticed that both CPI and Gaps-statistics does not decline too much at \(k\) of 5. Therefore, considering these knowledge, \(k\) of 5 is chosen as the optimal clustering number for futher analyses.

4) get results from single algorithm

In this part I will first show how to use MOVICS to perform multi-omics integrative clustering by specifying one algorithm with detailed parameters. For example, let us try iClusterBayes method like below:

Otherwise, a unified function can be used for all algorithms individually with detailed parameters, that is, getMOIC().

By specifying only one algorithm (i.e., iClusterBayes) in the argument of methodslist, the above getMOIC() will return exactly the same results from getiClusterBayes() if the same parameters are provided. The returned result contains a clust.res object that has two columns: clust to indicate the subtype that the sample belongs to, and samID records the corresponding sample name. For algorithms that provide feature selection procedure (i.e., iClusterBayes, CIMLR, and MoCluster), the result also contains a feat.res object that stores the information of such procedure. To those algorithms involving hierarchical clustering (e.g., COCA, ConsensusClustering), the corresponding dendrogram for sample clustering will be also returned as clust.dend.

5) get results from multiple algorithms at once

If you simultaneously specify a list of algorithms to methodslist argument in getMOIC(), it will automatically perform each algorithm with default parameters one by one, and a list of results derived from specified algorithms will be finally returned. Now that iClusterBayes has been finished, let us try other 9 algorithms at once with parameters by default. This will take some time, so have a coffee break.

6) get consensus from different algorithms

Since now all the clustering results from 10 algorthms are in hand, MOVICS borrows the idea of consensus clustering for later integration of the clustering results derived from different algorithms, so as to improve the clustering robustness. To be specific, if \(t_{max}\) algorithms are specified where \(2\le t_{max} \le10\), getConsensusMOIC() calculates a matrix \(M_{n\times n}^{(t)}\) per algorithm where \(n\) is the number of samples, and \(M_{ij}^{(t)}=1\) only when the sample \(i\) and \(j\) are clustered in the same subtype, otherwise \(M_{ij}^{(t)}=0\). After get all results from specified algorithms, MOVICS calculates a consensus matrix \(CM=\sum_{t=1}^{t_{max}}M^{(t)}\), and \(cm_{ij}\in[0,10]\). Applying hierarchical clustering to \(CM\), MOVICS generates a clustering result which represents a robust pair-wise similarities for samples because it considers different multi-omics integrative clustering algorithms. The simplest way to perform getConsensusMOIC() is to pass a list of object returned by multiple get%algorithm_name%() or by getMOIC() with specific argument of methodslist. This can be done like below:

Figure 2. Consensus heatmap based on results from 10 multi-omics integrative clustering algorithms with cluster number of 5.

Figure 2. Consensus heatmap based on results from 10 multi-omics integrative clustering algorithms with cluster number of 5.

Remember, you must choose at least two methods to perform the above consensus clustering by getConsensusMOIC(), and this function will return a consensus matrix (a probability matrix represents how many times samples belonging to the same subtype can be clustered together by different multi-omics clustering methods) and a corresponding consensus heatmap. Ideally, the consensus heatmap will show a perfect diagonal rectangle, and the input values are 0 and 1 only because all algorithms derived the same clustering results.

7) get multi-omics heatmap based on clustering result

Genome-wide heatmaps are widely used to graphically display potential underlying patterns within the large genomic dataset. They have been used to reveal information about how the samples/genes cluster together and provide insights into potential sample biases or other artifacts. Herein MOVICS provides getMoHeatmap to visually deal with preclustered multi-omics data and generate an exquisite heatmap for publication requirement. Before using getMoHeatmap(), omics data should be properly processed by using the function of getStdiz() which returns a list storing normalized omics data. Omics data, especially for expression (e.g., RNA and protein), should be centered (centerFlag = T) or scaled (scaleFlag = T) or z-scored (both centered and scaled). Generally, DNA methylation data (\(\beta\) matrix ranging from 0 to 1) and somatic mutation (0 and 1 binary matrix) should not be normalized. However, it is a good choice to convert the methylation \(\beta\) value to M value following the formula of \(M=log_2\frac{\beta}{1-\beta}\) for its stronger signal in visualization and M value is suitable for normalization. This function also provides an argument of halfwidth for continuous omics data; such argument is used to truncate the ‘extremum’ after normalization; specifically, normalized values that exceed the halfwidth boundaries will be replaced by the halfwidth, which is vary useful to map colors in heatmap.

As I mentioned earlier, several algorithms also provide feature selection; those selected features show a complex cross-talk with other omics data and might have special biological significance that drives the heterogeneity of cancers. Therefore I show below how to generate a comprehensive heatmap based on a single algorithm (e.g., iClusterBayes) with selected features by using getMoHeatmap(). However in the first place, features must be extracted:

The feat.res contained in iClusterBayes.res is sorted by posterior probability of features for each omics data. In this manner, the top 10 features are selected for each omics data and a feature list is generated and named as annRow for heatmap raw annotation (Do not worry, I will talk about how to attach column annotation for samples in a minute).

Figure 3. Comprehensive heatmap of multi-omics integrative clustering by iClusterBayes algorithm with annotation of potential cluster drivers (selected features).

Figure 3. Comprehensive heatmap of multi-omics integrative clustering by iClusterBayes algorithm with annotation of potential cluster drivers (selected features).

Of course, since there are a total of 10 results stored in the moic.res.list, you can also choose any one of them to create the heatmap. Here I select COCA which also returns sample dendrogram like below:

Figure 4. Comprehensive heatmap of multi-omics integrative clustering by COCA algorithm with dendrogram for samples.

Figure 4. Comprehensive heatmap of multi-omics integrative clustering by COCA algorithm with dendrogram for samples.

Now go back to the consensus result of cmoic.brca that integrates 10 algorithms and this time samples’ annotation is also provided to generate the heatmap. Since the core function of getMoHeatmap() is based on ComplexHeatmap R package, when creating annotations, you should always use circlize::colorRamp2() function to generate the color mapping function for continuous variables (e.g., age in this example).

Figure 5. Comprehensive heatmap of multi-omics integrative clustering based on consensus across 10 algorithms with clinicopathological annotation for samples.

Figure 5. Comprehensive heatmap of multi-omics integrative clustering based on consensus across 10 algorithms with clinicopathological annotation for samples.

4.2.2 COMP Module

After identification of cancer subtypes, it is essential to further characterize each subtype by discovering their difference from multile aspects. To this end, MOVICS provides commonly used downstream analyses in cancer subtyping researches for easily cohesion with results derived from GET Module. Now let us check them out.

1) compare survival outcome

I start with comparing the prognosis among different subtypes. MOVICS provides function of compSurv() which not only calculates the overall nominal P value by log-rank test, but also performs pair-wise comparsion and derives adjusted P values if more than two subtypes are identified. These information will be all printed in the Kaplan-Meier Curve which is convenient for researchers to refer. Except for clustering results (e.g., cmoic.brca in this example), you must additionally provide surv.info argument which should be a data.frame that stores a futime column for survival time (unit of day) and another fustat column for survival outcome (0 = censor; 1 = event).

Figure 6. Kaplan-Meier survival curve of 5 identified subtypes of breast cancer in TCGA-BRCA cohort.

Figure 6. Kaplan-Meier survival curve of 5 identified subtypes of breast cancer in TCGA-BRCA cohort.

2) compare clinical features

Then compare the clinical features among different subtypes. MOVICS provides function of compClinvar() which enables summarizing both continous and categorical variables and performing proper statistical tests. This function can give a table in .docx format that is easy to use in medical research papers.

Table 1. Comparison of clinical features among 5 identified subtype of breast cancer in TCGA-BRCA cohort.
level CS1 CS2 CS3 CS4 CS5 p test sig
n 146 133 107 144 113
fustat (%) 0 118 (80.8) 115 (86.5) 100 (93.5) 137 (95.1) 97 (85.8) 0.001 **
1 28 (19.2) 18 (13.5) 7 ( 6.5) 7 ( 4.9) 16 (14.2)
futime (median [IQR]) 719.00 [431.50, 1355.00] 746.50 [432.25, 1420.25] 703.00 [436.00, 1270.50] 848.50 [470.50, 1547.75] 742.00 [470.00, 1605.00] 0.691 nonnorm
PAM50 (%) Basal 2 ( 1.4) 0 ( 0.0) 0 ( 0.0) 0 ( 0.0) 109 (96.5) <0.001 **
Her2 38 (26.0) 0 ( 0.0) 0 ( 0.0) 0 ( 0.0) 0 ( 0.0)
LumA 47 (32.2) 78 (58.6) 91 (85.0) 128 (88.9) 0 ( 0.0)
LumB 53 (36.3) 55 (41.4) 16 (15.0) 0 ( 0.0) 0 ( 0.0)
Normal 6 ( 4.1) 0 ( 0.0) 0 ( 0.0) 16 (11.1) 4 ( 3.5)
pstage (%) T1 40 (27.4) 35 (26.3) 32 (29.9) 40 (27.8) 25 (22.1) NA exact
T2 85 (58.2) 77 (57.9) 63 (58.9) 72 (50.0) 70 (61.9)
T3 12 ( 8.2) 17 (12.8) 10 ( 9.3) 31 (21.5) 14 (12.4)
T4 9 ( 6.2) 4 ( 3.0) 2 ( 1.9) 1 ( 0.7) 3 ( 2.7)
TX 0 ( 0.0) 0 ( 0.0) 0 ( 0.0) 0 ( 0.0) 1 ( 0.9)
age 58.47 (13.75) 59.32 (14.03) 59.62 (12.40) 57.58 (12.23) 55.67 (12.11) 0.116

3) compare mutational frequency

Subype-specific mutation might be promising as therapeutic target. Hence, I then compare the mutational frequency among different clusters. MOVICS provides function of compMut() which deals with binary mutational data (0 = wild; 1 = mutated). This function applies independent testing (i.e., Fisher' s exact test or \(\chi^2\) test) for each mutation, and generates a table with statistical results (.docx file also if specified), and creates an oncoprint using those differentially mutated genes.

Figure 7. Mutational oncoprint of 5 identified subtypes of breast cancer in TCGA-BRCA cohort.

Figure 7. Mutational oncoprint of 5 identified subtypes of breast cancer in TCGA-BRCA cohort.

Table 2. Comparison of mutational frequency among 5 identified subtype of breast cancer in TCGA-BRCA cohort.
Gene (Mutated) TMB CS1 CS2 CS3 CS4 CS5 pvalue padj
PIK3CA 208 (32%) 55 (37.7%) 10 ( 7.5%) 77 (72.0%) 64 (44.4%) 2 ( 1.8%) 6.90e-42 3.10e-41
TP53 186 (29%) 83 (56.8%) 12 ( 9.0%) 0 ( 0.0%) 8 ( 5.6%) 83 (73.5%) 2.06e-63 1.85e-62
TTN 111 (17%) 41 (28.1%) 14 (10.5%) 14 (13.1%) 22 (15.3%) 20 (17.7%) 2.03e-03 3.04e-03
CDH1 83 (13%) 11 ( 7.5%) 3 ( 2.3%) 19 (17.8%) 49 (34.0%) 1 ( 0.9%) 5.38e-19 1.61e-18
GATA3 58 ( 9%) 5 ( 3.4%) 32 (24.1%) 10 ( 9.3%) 11 ( 7.6%) 0 ( 0.0%) 6.21e-11 1.40e-10
MLL3 49 ( 8%) 11 (7.5%) 12 (9.0%) 10 (9.3%) 11 (7.6%) 5 (4.4%) 6.31e-01 6.31e-01
MUC16 48 ( 8%) 18 (12.3%) 8 ( 6.0%) 6 ( 5.6%) 8 ( 5.6%) 8 ( 7.1%) 2.05e-01 2.31e-01
MAP3K1 38 ( 6%) 2 ( 1.4%) 6 ( 4.5%) 16 (15.0%) 12 ( 8.3%) 2 ( 1.8%) 3.65e-05 6.57e-05
SYNE1 33 ( 5%) 9 (6.2%) 1 (0.8%) 5 (4.7%) 8 (5.6%) 10 (8.8%) 3.27e-02 4.20e-02

4) compare total mutation burden

Needless to say, immunotherapy is becoming a pillar of modern cancer treatment. Recent analyses have linked the tumoral genomic landscape with antitumor immunity. In particular, it has been proposed that overall mutational load might drive T-cell responses7, whereas tumor aneuploidy correlates with markers of immune evasion and reduced response to immunotherapy8. To quantify these genomic alterations that may affect immunotherapy, MOVICS provides two functions to calculate total mutation burden (TMB) and fraction genome altered (FGA). To be specific, Total Mutations Burden (TMB) refers to the number of mutations that are found in the tumor genome, while FGA is the percentage of genome that has been affected by copy number gains or losses . Both attributes are useful for genetic researchers as they provide them with more in-depth information on the genomic make-up of the tumors. Let me show you how to use these two functions by starting with compTMB(). First of all, the input maf data for this function must have the following 10 columns at least:

Run compTMB() then. By default, this function considers nonsynonymous variants only when counting somatic mutation frequency, including Frame_Shift_Del, Frame_Shift_Ins, Splice_Site, Translation_Start_Site, Nonsense_Mutation, Nonstop_Mutation, In_Frame_Del, In_Frame_Ins, and Missense_Mutation by default. In addition to TMB, this function also classifies Single Nucleotide Variants into Transitions and Transversions (TiTv), and depicts the distribution of both TMB and TiTv.

Figure 8. Comparison of TMB and TiTv among 5 identified subtypes of breast cancer in TCGA-BRCA cohort.

Figure 8. Comparison of TMB and TiTv among 5 identified subtypes of breast cancer in TCGA-BRCA cohort.

Table 3. Demo of comparison of TMB among 5 identified subtype of breast cancer in TCGA-BRCA cohort.

5) compare fraction genome altered

Next, compFGA() calculates not only FGA but also computes specific gain (FGG) or loss (FGL) per sample within each subtype. To get this function worked, an eligible input of segmented copy number should be prepared with exactly the same column name like below:

Let’s see how the input looks like:

Then this function can be run without any difficulties. Notably, if your CNA calling procedure did not provide a segmented copy number as value column but the original copy number, argument of iscopynumber must be switched to TRUE instead.

Figure 9. Barplot of fraction genome altered among 5 identified subtypes of breast cancer in TCGA-BRCA cohort.

Figure 9. Barplot of fraction genome altered among 5 identified subtypes of breast cancer in TCGA-BRCA cohort.

Table 4. Demo of comparison of fraction genome altered among 5 identified subtype of breast cancer in TCGA-BRCA cohort.
samID FGA FGG FGL Subtype
BRCA-A03L-01A 0.6217991 0.3086727 0.3131264 CS1
BRCA-A04R-01A 0.2531019 0.0913201 0.1617818 CS2
BRCA-A075-01A 0.7007067 0.4144424 0.2862643 CS1
BRCA-A08O-01A 0.6501287 0.4564814 0.1936473 CS3
BRCA-A0A6-01A 0.1468893 0.0635649 0.0833244 CS4
BRCA-A0AD-01A 0.1722214 0.0386452 0.1335762 CS2

6) compare drug sensitivity

Predicting response to medication is particularly important for drugs with a narrow therapeutic index, for example chemotherapeutic agents, because response is highly variable and side effects are potentially lethal. Therefore, Paul Geeleher et al. (2014)9 used baseline gene expression and \(in\;vitro\) drug sensitivity derived from cell lines, coupled with \(in\;vivo\) baseline tumor gene expression, to predict patients’ response to drugs. Paul developped an R package pRRophetic for prediction of clinical chemotherapeutic response from tumor gene expression levels10, and now this function has been involved in MOVICS to examine difference of drug sensitivity among different clusters. Here I estimate the \(IC_{50}\) of Cisplatin for 5 identified subtypes of breast cancer in TCGA.

Figure 10. Boxviolins for estimated IC50 of Cisplatin and Paclitaxel among 5 identified subtypes of breast cancer in TCGA-BRCA cohort.Figure 10. Boxviolins for estimated IC50 of Cisplatin and Paclitaxel among 5 identified subtypes of breast cancer in TCGA-BRCA cohort.

Figure 10. Boxviolins for estimated IC50 of Cisplatin and Paclitaxel among 5 identified subtypes of breast cancer in TCGA-BRCA cohort.

Table 5. Demo of estimated IC50 for Cisplatin among 5 identified subtype of breast cancer in TCGA-BRCA cohort.
Est.IC50 Subtype
BRCA-A03L-01A 3.801326 CS1
BRCA-A075-01A 3.834292 CS1
BRCA-A0AW-01A 3.127552 CS1
BRCA-A0AZ-01A 3.302353 CS1
BRCA-A0BC-01A 3.409796 CS1
BRCA-A0BF-01A 3.357585 CS1

7) compare agreement with other subtypes

At present, many cancers have traditional classifications, and evaluating the consistency of new subtypes with previous classifications is critical to reflect the robustness of clustering analysis and to determine potential but novel subtypes. To measure the agreement (similarity) between the current subtypes and other pre-existed classifications MOVICS provides function of compAgree() to calculate four statistics: Rand Index (RI)11, Adjusted Mutual Information (AMI)12, Jaccard Index (JI)13, and Fowlkes-Mallows (FM)14; all these measurements range from 0 to 1 and the larger the value is, the more similar the two appraises are. This function can also generate an alluvial diagram to visualize the agreement of two appraises with the current subtypes as reference.

Figure 11. Agreement of 5 identified subtypes of breast cancer with PAM50 classification and pathological stage in TCGA-BRCA cohort.

Figure 11. Agreement of 5 identified subtypes of breast cancer with PAM50 classification and pathological stage in TCGA-BRCA cohort.

Table 6. Agreement of 5 identified subtypes with PAM50 classification and pathological stage in TCGA-BRCA cohort.
current.subtype other.subtype RI AMI JI FM
Subtype PAM50 0.6929744 0.4043456 0.2910887 0.4694389
Subtype pstage 0.5506751 0.0124902 0.1565996 0.2884961

4.2.3 RUN Module

Take a deep breath, we successfully enter the last RUN Module, and victory is in sight. In this module, MOVICS aims to characterize different subtypes by identifying their potential predictive biomarkers and functional pathways. Identifying and applying molecular biomarkers to predict subtype with efficiency is particularly important for disease management and treatment, thus improving clinical outcome. In this context, MOVICS searches for subtype-specific biomarkers by starting with differential expression analysis (DEA).

1) run differential expression analysis

MOVICS provides runDEA() function which embeds three state-of-the-art DEA approches to identify differentially expressed genes (DEGs), including edgeR and DESeq2 for RNA-Seq count data and limma for microarray data. Notably, since runDEA() checks the data scale automatically when choosing limma algorithm, it is recommended to provide a normalized expression data (e.g., RSEM, FPKM, TPM) witout z-score or log2 transformation. This step is also quite time-consuming, especially in the case of using raw count data, so release your hand from the keyboard, close your eyes, and relax for a moment.

Each identified cancer subtype will be compared with the rest (Others) and the corresponding .txt file will be stored according to argument of res.path. By default, these files will be saved under the current working directory.

2) run biomarker identification procedure

In this procedure, the most differentially expressed genes sorted by log2FoldChange are chosen as the biomarkers for each subtype. These biomarkers should pass the significance threshold (e.g., nominal \(P\) value < 0.05 and adjusted \(P\) value < 0.05) and must not overlap with any biomarkers indentified for other subtypes.

Figure 12. Heatmap of subtype-specific upregulated biomarkers using edgeR algorithm for 5 identified subtypes of breast cancer in TCGA-BRCA cohort.

Figure 12. Heatmap of subtype-specific upregulated biomarkers using edgeR algorithm for 5 identified subtypes of breast cancer in TCGA-BRCA cohort.

Table 7. Demo of subtype-specific upregulated biomarkers for 5 identified subtypes of breast cancer in TCGA-BRCA cohort.
probe class dirct
PNMT CS1 up
AKR1B15 CS1 up
DLK1 CS1 up
ACE2 CS1 up
CRISP3 CS1 up
ACSM1 CS1 up

Then try results derived from limma to identify subtype-specific downregulated biomarkers. Please feel free to try DESeq2.

Figure 13. Heatmap of subtype-specific downregulated biomarkers using limma algorithm for 5 identified subtypes of breast cancer in TCGA-BRCA cohort.

Figure 13. Heatmap of subtype-specific downregulated biomarkers using limma algorithm for 5 identified subtypes of breast cancer in TCGA-BRCA cohort.

3) run gene set enrichment analysis

Similarly, GSEA is run for each subtype based on its corresponding DEA result to identify subtype-specific functional pathways. To this end, I have prepared a geneset background which includes all gene sets derived from GO biological processes (c5.bp.v7.1.symbols.gmt) from The Molecular Signatures Database (MSigDB, https://www.gsea-msigdb.org/gsea/msigdb/index.jsp).

Likewise, these identified specific pathways should pass the significance threshold (e.g., nominal \(P\) value < 0.05 and adjusted \(P\) value < 0.25) and must not overlap with any pathways identified for other subtypes. After having the subtype-specific pathways, genes that are inside the pathways are retrieved to calculate a single sample enrichment score by using GSVA R package. Subsequently, subtype-specific enrichment score will be represented by the mean or median (mean by default) value within the subtype, and will be further visualized by diagonal heatmap.

# run GSEA to identify up-regulated GO pathways using results from edgeR
gsea.up <- runGSEA(moic.res     = cmoic.brca,
                   dea.method   = "edger", # name of DEA method
                   prefix       = "TCGA-BRCA", # MUST be the same of argument in runDEA()
                   dat.path     = getwd(), # path of DEA files
                   res.path     = getwd(), # path to save GSEA files
                   msigdb.path  = MSIGDB.FILE, # MUST be the ABSOLUTE path of msigdb file
                   norm.expr    = fpkm, # use normalized expression to calculate enrichment score
                   dirct        = "up", # direction of dysregulation in pathway
                   p.cutoff     = 0.05, # p cutoff to identify significant pathways
                   p.adj.cutoff = 0.25, # padj cutoff to identify significant pathways
                   gsva.method  = "gsva", # method to calculate single sample enrichment score
                   norm.method  = "mean", # normalization method to calculate subtype-specific enrichment score
                   fig.name     = "UPREGULATED PATHWAY HEATMAP")
#> --all samples matched.
#> Loading required package: GSEABase
#> Loading required package: annotate
#> Loading required package: AnnotationDbi
#> Loading required package: stats4
#> Loading required package: IRanges
#> Loading required package: S4Vectors
#> 
#> Attaching package: 'S4Vectors'
#> The following objects are masked from 'package:dplyr':
#> 
#>     first, rename
#> The following object is masked from 'package:base':
#> 
#>     expand.grid
#> 
#> Attaching package: 'IRanges'
#> The following objects are masked from 'package:dplyr':
#> 
#>     collapse, desc, slice
#> The following object is masked from 'package:grDevices':
#> 
#>     windows
#> 
#> Attaching package: 'AnnotationDbi'
#> The following object is masked from 'package:dplyr':
#> 
#>     select
#> Loading required package: XML
#> Loading required package: graph
#> 
#> Attaching package: 'graph'
#> The following object is masked from 'package:XML':
#> 
#>     addNode
#> GSEA done...
#> --log2 transformation done for expression data.
#> Estimating GSVA scores for 50 gene sets.
#> Computing observed enrichment scores
#> Estimating ECDFs with Gaussian kernels
#> Using parallel with 1 cores
#> 
  |                                                                       
  |                                                                 |   0%
  |                                                                       
  |=                                                                |   2%
  |                                                                       
  |===                                                              |   4%
  |                                                                       
  |====                                                             |   6%
  |                                                                       
  |=====                                                            |   8%
  |                                                                       
  |======                                                           |  10%
  |                                                                       
  |========                                                         |  12%
  |                                                                       
  |=========                                                        |  14%
  |                                                                       
  |==========                                                       |  16%
  |                                                                       
  |============                                                     |  18%
  |                                                                       
  |=============                                                    |  20%
  |                                                                       
  |==============                                                   |  22%
  |                                                                       
  |================                                                 |  24%
  |                                                                       
  |=================                                                |  26%
  |                                                                       
  |==================                                               |  28%
  |                                                                       
  |====================                                             |  30%
  |                                                                       
  |=====================                                            |  32%
  |                                                                       
  |======================                                           |  34%
  |                                                                       
  |=======================                                          |  36%
  |                                                                       
  |=========================                                        |  38%
  |                                                                       
  |==========================                                       |  40%
  |                                                                       
  |===========================                                      |  42%
  |                                                                       
  |=============================                                    |  44%
  |                                                                       
  |==============================                                   |  46%
  |                                                                       
  |===============================                                  |  48%
  |                                                                       
  |================================                                 |  50%
  |                                                                       
  |==================================                               |  52%
  |                                                                       
  |===================================                              |  54%
  |                                                                       
  |====================================                             |  56%
  |                                                                       
  |======================================                           |  58%
  |                                                                       
  |=======================================                          |  60%
  |                                                                       
  |========================================                         |  62%
  |                                                                       
  |==========================================                       |  64%
  |                                                                       
  |===========================================                      |  66%
  |                                                                       
  |============================================                     |  68%
  |                                                                       
  |==============================================                   |  70%
  |                                                                       
  |===============================================                  |  72%
  |                                                                       
  |================================================                 |  74%
  |                                                                       
  |=================================================                |  76%
  |                                                                       
  |===================================================              |  78%
  |                                                                       
  |====================================================             |  80%
  |                                                                       
  |=====================================================            |  82%
  |                                                                       
  |=======================================================          |  84%
  |                                                                       
  |========================================================         |  86%
  |                                                                       
  |=========================================================        |  88%
  |                                                                       
  |==========================================================       |  90%
  |                                                                       
  |============================================================     |  92%
  |                                                                       
  |=============================================================    |  94%
  |                                                                       
  |==============================================================   |  96%
  |                                                                       
  |================================================================ |  98%
  |                                                                       
  |=================================================================| 100%
#> gsva done...
#> heatmap done...
Figure 14. Heatmap of subtype-specific upregulated pathways using edgeR algorithm for 5 identified subtypes in TCGA-BRCA cohort.

Figure 14. Heatmap of subtype-specific upregulated pathways using edgeR algorithm for 5 identified subtypes in TCGA-BRCA cohort.

Check some columns of GSEA results for the first subtype (CS1).

Table 8. Demo of GSEA results for the first cancer subtype (CS1) of breast cancer in TCGA-BRCA cohort.
setSize enrichmentScore NES pvalue
GO_REGULATION_OF_ION_TRANSPORT 420 -0.4832370 -1.555051 0.0010040
GO_EMBRYONIC_MORPHOGENESIS 459 -0.4504464 -1.449967 0.0010050
GO_EPITHELIAL_CELL_DIFFERENTIATION 470 -0.4956895 -1.597430 0.0010050
GO_SYNAPTIC_SIGNALING 441 -0.5216787 -1.677297 0.0010050
GO_BEHAVIOR 392 -0.5060951 -1.625033 0.0010091
GO_REGULATION_OF_TRANS_SYNAPTIC_SIGNALING 293 -0.5122532 -1.631490 0.0010101

Also check results of subtype-specific enrichment scores.

Table 9. Demo of subtype-specific enrichment scores among 5 identified subtypes of breast cancer in TCGA-BRCA cohort.
CS1 CS2 CS3 CS4 CS5
GO_BENZENE_CONTAINING_COMPOUND_METABOLIC_PROCESS 0 0 0 0 -1
GO_REGULATION_OF_MONOCYTE_CHEMOTAXIS 0 0 0 0 0

Then try results derived from DESeq2 to identify subtype-specific downregulated pathways. Please feel free to try limma.

# run GSEA to identify down-regulated GO pathways using results from DESeq2
gsea.dn <- runGSEA(moic.res     = cmoic.brca,
                   dea.method   = "deseq2",
                   prefix       = "TCGA-BRCA",
                   msigdb.path  = MSIGDB.FILE,
                   norm.expr    = fpkm,
                   dirct        = "down",
                   p.cutoff     = 0.05,
                   p.adj.cutoff = 0.25,
                   gsva.method  = "ssgsea", # switch to ssgsea
                   norm.method  = "median", # switch to median
                   fig.name     = "DOWNREGULATED PATHWAY HEATMAP") 
#> --all samples matched.
#> GSEA done...
#> --log2 transformation done for expression data.
#> Estimating ssGSEA scores for 50 gene sets.
#> 
  |                                                                       
  |                                                                 |   0%Using parallel with 1 cores
#> 
  |                                                                       
  |                                                                 |   1%
  |                                                                       
  |=                                                                |   1%
  |                                                                       
  |=                                                                |   2%
  |                                                                       
  |==                                                               |   2%
  |                                                                       
  |==                                                               |   3%
  |                                                                       
  |==                                                               |   4%
  |                                                                       
  |===                                                              |   4%
  |                                                                       
  |===                                                              |   5%
  |                                                                       
  |====                                                             |   5%
  |                                                                       
  |====                                                             |   6%
  |                                                                       
  |====                                                             |   7%
  |                                                                       
  |=====                                                            |   7%
  |                                                                       
  |=====                                                            |   8%
  |                                                                       
  |======                                                           |   9%
  |                                                                       
  |======                                                           |  10%
  |                                                                       
  |=======                                                          |  10%
  |                                                                       
  |=======                                                          |  11%
  |                                                                       
  |=======                                                          |  12%
  |                                                                       
  |========                                                         |  12%
  |                                                                       
  |========                                                         |  13%
  |                                                                       
  |=========                                                        |  13%
  |                                                                       
  |=========                                                        |  14%
  |                                                                       
  |==========                                                       |  15%
  |                                                                       
  |==========                                                       |  16%
  |                                                                       
  |===========                                                      |  16%
  |                                                                       
  |===========                                                      |  17%
  |                                                                       
  |===========                                                      |  18%
  |                                                                       
  |============                                                     |  18%
  |                                                                       
  |============                                                     |  19%
  |                                                                       
  |=============                                                    |  19%
  |                                                                       
  |=============                                                    |  20%
  |                                                                       
  |=============                                                    |  21%
  |                                                                       
  |==============                                                   |  21%
  |                                                                       
  |==============                                                   |  22%
  |                                                                       
  |===============                                                  |  22%
  |                                                                       
  |===============                                                  |  23%
  |                                                                       
  |===============                                                  |  24%
  |                                                                       
  |================                                                 |  24%
  |                                                                       
  |================                                                 |  25%
  |                                                                       
  |=================                                                |  26%
  |                                                                       
  |=================                                                |  27%
  |                                                                       
  |==================                                               |  27%
  |                                                                       
  |==================                                               |  28%
  |                                                                       
  |===================                                              |  29%
  |                                                                       
  |===================                                              |  30%
  |                                                                       
  |====================                                             |  30%
  |                                                                       
  |====================                                             |  31%
  |                                                                       
  |=====================                                            |  32%
  |                                                                       
  |=====================                                            |  33%
  |                                                                       
  |======================                                           |  33%
  |                                                                       
  |======================                                           |  34%
  |                                                                       
  |======================                                           |  35%
  |                                                                       
  |=======================                                          |  35%
  |                                                                       
  |=======================                                          |  36%
  |                                                                       
  |========================                                         |  36%
  |                                                                       
  |========================                                         |  37%
  |                                                                       
  |========================                                         |  38%
  |                                                                       
  |=========================                                        |  38%
  |                                                                       
  |=========================                                        |  39%
  |                                                                       
  |==========================                                       |  39%
  |                                                                       
  |==========================                                       |  40%
  |                                                                       
  |==========================                                       |  41%
  |                                                                       
  |===========================                                      |  41%
  |                                                                       
  |===========================                                      |  42%
  |                                                                       
  |============================                                     |  42%
  |                                                                       
  |============================                                     |  43%
  |                                                                       
  |============================                                     |  44%
  |                                                                       
  |=============================                                    |  44%
  |                                                                       
  |=============================                                    |  45%
  |                                                                       
  |==============================                                   |  45%
  |                                                                       
  |==============================                                   |  46%
  |                                                                       
  |==============================                                   |  47%
  |                                                                       
  |===============================                                  |  47%
  |                                                                       
  |===============================                                  |  48%
  |                                                                       
  |================================                                 |  49%
  |                                                                       
  |================================                                 |  50%
  |                                                                       
  |=================================                                |  50%
  |                                                                       
  |=================================                                |  51%
  |                                                                       
  |==================================                               |  52%
  |                                                                       
  |==================================                               |  53%
  |                                                                       
  |===================================                              |  53%
  |                                                                       
  |===================================                              |  54%
  |                                                                       
  |===================================                              |  55%
  |                                                                       
  |====================================                             |  55%
  |                                                                       
  |====================================                             |  56%
  |                                                                       
  |=====================================                            |  56%
  |                                                                       
  |=====================================                            |  57%
  |                                                                       
  |=====================================                            |  58%
  |                                                                       
  |======================================                           |  58%
  |                                                                       
  |======================================                           |  59%
  |                                                                       
  |=======================================                          |  59%
  |                                                                       
  |=======================================                          |  60%
  |                                                                       
  |=======================================                          |  61%
  |                                                                       
  |========================================                         |  61%
  |                                                                       
  |========================================                         |  62%
  |                                                                       
  |=========================================                        |  62%
  |                                                                       
  |=========================================                        |  63%
  |                                                                       
  |=========================================                        |  64%
  |                                                                       
  |==========================================                       |  64%
  |                                                                       
  |==========================================                       |  65%
  |                                                                       
  |===========================================                      |  65%
  |                                                                       
  |===========================================                      |  66%
  |                                                                       
  |===========================================                      |  67%
  |                                                                       
  |============================================                     |  67%
  |                                                                       
  |============================================                     |  68%
  |                                                                       
  |=============================================                    |  69%
  |                                                                       
  |=============================================                    |  70%
  |                                                                       
  |==============================================                   |  70%
  |                                                                       
  |==============================================                   |  71%
  |                                                                       
  |===============================================                  |  72%
  |                                                                       
  |===============================================                  |  73%
  |                                                                       
  |================================================                 |  73%
  |                                                                       
  |================================================                 |  74%
  |                                                                       
  |=================================================                |  75%
  |                                                                       
  |=================================================                |  76%
  |                                                                       
  |==================================================               |  76%
  |                                                                       
  |==================================================               |  77%
  |                                                                       
  |==================================================               |  78%
  |                                                                       
  |===================================================              |  78%
  |                                                                       
  |===================================================              |  79%
  |                                                                       
  |====================================================             |  79%
  |                                                                       
  |====================================================             |  80%
  |                                                                       
  |====================================================             |  81%
  |                                                                       
  |=====================================================            |  81%
  |                                                                       
  |=====================================================            |  82%
  |                                                                       
  |======================================================           |  82%
  |                                                                       
  |======================================================           |  83%
  |                                                                       
  |======================================================           |  84%
  |                                                                       
  |=======================================================          |  84%
  |                                                                       
  |=======================================================          |  85%
  |                                                                       
  |========================================================         |  86%
  |                                                                       
  |========================================================         |  87%
  |                                                                       
  |=========================================================        |  87%
  |                                                                       
  |=========================================================        |  88%
  |                                                                       
  |==========================================================       |  88%
  |                                                                       
  |==========================================================       |  89%
  |                                                                       
  |==========================================================       |  90%
  |                                                                       
  |===========================================================      |  90%
  |                                                                       
  |===========================================================      |  91%
  |                                                                       
  |============================================================     |  92%
  |                                                                       
  |============================================================     |  93%
  |                                                                       
  |=============================================================    |  93%
  |                                                                       
  |=============================================================    |  94%
  |                                                                       
  |=============================================================    |  95%
  |                                                                       
  |==============================================================   |  95%
  |                                                                       
  |==============================================================   |  96%
  |                                                                       
  |===============================================================  |  96%
  |                                                                       
  |===============================================================  |  97%
  |                                                                       
  |===============================================================  |  98%
  |                                                                       
  |================================================================ |  98%
  |                                                                       
  |================================================================ |  99%
  |                                                                       
  |=================================================================|  99%
  |                                                                       
  |=================================================================| 100%
#> ssgsea done...
#> heatmap done...
Figure 15. Heatmap of subtype-specific downregulated pathways using limma algorithm for 5 identified subtypes in TCGA-BRCA cohort.

Figure 15. Heatmap of subtype-specific downregulated pathways using limma algorithm for 5 identified subtypes in TCGA-BRCA cohort.

4) run nearest template prediction in external cohort

Oh wait, did we forget something? Yeah there is a dataset that has not been used yet, so let’s see if we can use those subtype-specific biomarkers to verify the current breast cancer subtypes in the external Yau cohort. In this part, our core purpose is to predict the possible subtypes of each sample in the external dataset. In most cases, this is a multi-classification problem, and the identified biomarkers may be difficult to be overall matched in the external cohort, so it might not be reliable to use model-based predictive algorithms. Therefore, MOVICS switches to nearest template prediction (NTP) which can be flexibly applied to cross-platform, cross-species, and multiclass predictions without any optimization of analysis parameters15. Only one thing have to do is generating a template file, which has been fortunately prepared already.

Figure 16. Heatmap of NTP in Yau cohort using subtype-specific upregulated biomarkers identified from TCGA-BRCA cohort

Figure 16. Heatmap of NTP in Yau cohort using subtype-specific upregulated biomarkers identified from TCGA-BRCA cohort

#> 
#>  CS1  CS2  CS3  CS4  CS5 <NA> 
#>  104   85   90  120  140  143
Table 10. Demo of predicted subtypes in Yau cohort by NTP using subtype-specific upregulated biomarkers identified from TCGA-BRCA cohort.
prediction d.CS1 d.CS2 d.CS3 d.CS4 d.CS5 p.value FDR
107 CS2 0.7344 0.5165 0.7377 0.7471 0.7512 0.001 0.0020
109 CS1 0.5816 0.7581 0.7172 0.7308 0.7295 0.001 0.0020
11 CS1 0.6527 0.7307 0.7430 0.7557 0.7731 0.001 0.0020
110 CS2 0.7619 0.5418 0.7583 0.7655 0.7567 0.001 0.0020
111 CS1 0.6806 0.7157 0.7198 0.7889 0.7911 0.007 0.0106
113 CS4 0.6785 0.7505 0.7394 0.6389 0.7308 0.007 0.0106

Above an object brca.pred is prepared which is similar in structure of object returned by getMOIC() but only stores clust.res which can be passed to functions within COMP Module if additional data is available. For example, I herein first compare the survival outcome of the predicted 5 cancer subtypes in Yau cohort.

Figure 17. Kaplan-Meier survival curve of predicted 5 subtypes of breast cancer in Yau cohort.

Figure 17. Kaplan-Meier survival curve of predicted 5 subtypes of breast cancer in Yau cohort.

I further check the agreement between the predicted subtype and PAM50 classification.

Figure 18. Agreement of predicted 5 subtypes of breast cancer with PAM50 classification in Yau cohort.

Figure 18. Agreement of predicted 5 subtypes of breast cancer with PAM50 classification in Yau cohort.

Table 11. Agreement of 5 predicted subtypes of breast cancer with PAM50 classification in Yau cohort.
current.subtype other.subtype RI AMI JI FM
Subtype PAM50 0.7830213 0.3864876 0.3318526 0.4994425

It is clearly that the predicted subtypes of breast cancer in Yau cohort can distinguish prognosis well and also show similar pattern of agreement with PAM50 classification as compared to TCGA-BRCA to some extent.

5. Little Trick

Indeed, the moic.res parameter is necessary for almost all the functions, especially for downstream analyses. You may deem that if the moic.res object is not obtained properly through the module of GET, downstream analyses cannot run smoothly. But in fact, you can fool downstream modules (COMP and RUN) with little tricks. The core information stored in moic.res is a data.frame named clust.res with a samID column (character) for samples name (should be exactly the same with rownames) and a clust column (integer) for subtype indicator. Another information requried for some functions in downstream analyses is mo.method, a stringe value that is usually considered as prefix for output files. For example, if I want to check the prognostic value of PAM50 classification, the only thing I have to prepare is a pseudo moic.res object with a mo.method just like below:

Let check how the pseudo moic.res object looks like.

Ok everything is ready, just keep in mind how your interested subtypes are mapped and feel free to use such pseudo object to perform downstream analyses such as compSurv() as follows:

Figure 19. Kaplan-Meier survival curve of PAM50 subtypes of breast cancer with pseudo input in TCGA-BRCA cohort.

Figure 19. Kaplan-Meier survival curve of PAM50 subtypes of breast cancer with pseudo input in TCGA-BRCA cohort.

6. Summary

I have introduced the pipeline and described the functions included in MOVICS in almost detail. This package is rather young and I hope to improve it further in the future, along with more functionality. If you have any questions, bug reports, or suggestions for improving MOVICS, please email them to xlu.cpu@foxmail.com. Feeback is crucial for improving a package that tries to seamlessly incorporate functionality and flexibility from many other useful tools.

7. Session Information

sessionInfo()
#> R version 3.6.1 (2019-07-05)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 10 x64 (build 18362)
#> 
#> Matrix products: default
#> 
#> locale:
#> [1] LC_COLLATE=C                              
#> [2] LC_CTYPE=Chinese (Simplified)_China.936   
#> [3] LC_MONETARY=Chinese (Simplified)_China.936
#> [4] LC_NUMERIC=C                              
#> [5] LC_TIME=Chinese (Simplified)_China.936    
#> 
#> attached base packages:
#> [1] stats4    parallel  stats     graphics  grDevices utils     datasets 
#> [8] methods   base     
#> 
#> other attached packages:
#>  [1] GSEABase_1.48.0      graph_1.64.0         annotate_1.64.0     
#>  [4] XML_3.98-1.20        AnnotationDbi_1.48.0 IRanges_2.20.0      
#>  [7] S4Vectors_0.24.0     MOVICS_1.0.0         testthat_2.2.1      
#> [10] Biobase_2.46.0       BiocGenerics_0.32.0  devtools_2.2.1      
#> [13] usethis_1.5.1        kableExtra_1.1.2     knitr_1.25          
#> [16] dplyr_0.8.3         
#> 
#> loaded via a namespace (and not attached):
#>   [1] Hmisc_4.2-0                 corpcor_1.6.9              
#>   [3] aricode_1.0.0               class_7.3-15               
#>   [5] ps_1.3.0                    foreach_1.4.7              
#>   [7] rprojroot_1.3-2             crayon_1.3.4               
#>   [9] MASS_7.3-51.6               nlme_3.1-140               
#>  [11] backports_1.1.5             sva_3.34.0                 
#>  [13] impute_1.60.0               GOSemSim_2.12.0            
#>  [15] rlang_0.4.7                 readxl_1.3.1               
#>  [17] svd_0.5                     XVector_0.26.0             
#>  [19] heatmap.plus_1.3            irlba_2.3.3                
#>  [21] nloptr_1.2.1                callr_3.3.2                
#>  [23] limma_3.42.0                BiocParallel_1.20.0        
#>  [25] rjson_0.2.20                bit64_0.9-7                
#>  [27] glue_1.3.1                  rngtools_1.4               
#>  [29] SNFtool_2.3.0               processx_3.4.1             
#>  [31] NEMO_0.1.0                  oompaData_3.1.1            
#>  [33] DOSE_3.12.0                 haven_2.3.1                
#>  [35] tidyselect_0.2.5            SummarizedExperiment_1.16.0
#>  [37] km.ci_0.5-2                 rio_0.5.16                 
#>  [39] tidyr_1.0.0                 zoo_1.8-8                  
#>  [41] ggpubr_0.2.3                maftools_2.2.0             
#>  [43] ridge_2.4                   xtable_1.8-4               
#>  [45] magrittr_1.5                prettyGraphs_2.1.6         
#>  [47] evaluate_0.14               bibtex_0.4.2               
#>  [49] ggplot2_3.3.2               cli_1.1.0                  
#>  [51] zlibbioc_1.32.0             rstudioapi_0.10            
#>  [53] rpart_4.1-15                wordcloud_2.6              
#>  [55] fastmatch_1.1-0             ClassDiscovery_3.3.12      
#>  [57] InterSIM_2.2.0              shiny_1.4.0                
#>  [59] GSVA_1.34.0                 prettydoc_0.4.0            
#>  [61] xfun_0.10                   coca_1.1.0                 
#>  [63] clue_0.3-57                 pkgbuild_1.0.6             
#>  [65] cluster_2.1.0               caTools_1.17.1.2           
#>  [67] tidygraph_1.1.2             tibble_3.0.3               
#>  [69] flexclust_1.4-0             ggrepel_0.8.1              
#>  [71] png_0.1-7                   permute_0.9-5              
#>  [73] withr_2.1.2                 bitops_1.0-6               
#>  [75] ggforce_0.3.1               cellranger_1.1.0           
#>  [77] plyr_1.8.6                  e1071_1.7-2                
#>  [79] survey_4.0                  pillar_1.4.6               
#>  [81] RcppParallel_4.4.4          geepack_1.3-1              
#>  [83] gplots_3.0.1.1              GlobalOptions_0.1.1        
#>  [85] fs_1.3.1                    iClusterPlus_1.22.0        
#>  [87] GetoptLong_0.1.7            graphite_1.32.0            
#>  [89] europepmc_0.3               clusterProfiler_3.14.0     
#>  [91] vctrs_0.3.2                 ellipsis_0.3.0             
#>  [93] generics_0.0.2              urltools_1.7.3             
#>  [95] NMF_0.21.0                  tools_3.6.1                
#>  [97] foreign_0.8-71              entropy_1.2.1              
#>  [99] munsell_0.5.0               tweenr_1.0.1               
#> [101] fgsea_1.12.0                DelayedArray_0.12.0        
#> [103] abind_1.4-5                 fastmap_1.0.1              
#> [105] compiler_3.6.1              pkgload_1.0.2              
#> [107] httpuv_1.5.2                sessioninfo_1.1.1          
#> [109] pkgmaker_0.27               GenomeInfoDbData_1.2.2     
#> [111] gridExtra_2.3               edgeR_3.28.0               
#> [113] lattice_0.20-38             jstable_0.9.5              
#> [115] later_1.0.0                 alluvial_0.1-2             
#> [117] jsonlite_1.6                scales_1.1.1               
#> [119] dendsort_0.3.3              carData_3.0-2              
#> [121] pbapply_1.4-2               genefilter_1.68.0          
#> [123] promises_1.1.0              car_3.0-4                  
#> [125] SimDesign_2.0.1             doParallel_1.0.15          
#> [127] latticeExtra_0.6-28         checkmate_1.9.4            
#> [129] openxlsx_4.1.2              rmarkdown_2.3              
#> [131] cowplot_1.0.0               statmod_1.4.32             
#> [133] webshot_0.5.1               forcats_0.4.0              
#> [135] igraph_1.2.4.1              survival_3.2-3             
#> [137] IntNMF_1.2.0                yaml_2.2.0                 
#> [139] htmltools_0.4.0             memoise_1.1.0              
#> [141] modeltools_0.2-22           locfit_1.5-9.1             
#> [143] graphlayouts_0.5.0          viridisLite_0.3.0          
#> [145] digest_0.6.22               assertthat_0.2.1           
#> [147] mime_0.7                    rappdirs_0.3.1             
#> [149] registry_0.5-1              LRAcluster_1.0             
#> [151] KMsurv_0.1-5                RSQLite_2.1.2              
#> [153] remotes_2.1.0               data.table_1.12.6          
#> [155] blob_1.2.0                  vegan_2.5-6                
#> [157] mogsa_1.20.0                CMScaller_0.99.1           
#> [159] preprocessCore_1.48.0       survMisc_0.5.5             
#> [161] labeling_0.3                shinythemes_1.1.2          
#> [163] splines_3.6.1               Formula_1.2-3              
#> [165] RCurl_1.95-4.12             broom_0.7.0                
#> [167] hms_0.5.2                   colorspace_1.4-1           
#> [169] ConsensusClusterPlus_1.50.0 base64enc_0.1-3            
#> [171] BiocManager_1.30.9          GenomicRanges_1.38.0       
#> [173] shape_1.4.4                 nnet_7.3-12                
#> [175] Rcpp_1.0.2                  mclust_5.4.5               
#> [177] circlize_0.4.8              enrichplot_1.6.0           
#> [179] R6_2.4.0                    survminer_0.4.6            
#> [181] grid_3.6.1                  ggridges_0.5.1             
#> [183] lifecycle_0.2.0             acepack_1.4.1              
#> [185] labelled_2.5.0              oompaBase_3.2.9            
#> [187] zip_2.0.4                   curl_4.2                   
#> [189] ggsignif_0.6.0              minqa_1.2.4                
#> [191] gdata_2.18.0                PINSPlus_2.0.4             
#> [193] DO.db_2.9                   Matrix_1.2-17              
#> [195] qvalue_2.18.0               desc_1.2.0                 
#> [197] RColorBrewer_1.1-2          iterators_1.0.12           
#> [199] stringr_1.4.0               htmlwidgets_1.5.1          
#> [201] officer_0.3.12              polyclip_1.10-0            
#> [203] triebeard_0.3.0             purrr_0.3.3                
#> [205] CIMLR_1.0.0                 gridGraphics_0.4-1         
#> [207] rvest_0.3.4                 ComplexHeatmap_2.5.4       
#> [209] mgcv_1.8-28                 patchwork_0.0.1            
#> [211] htmlTable_1.13.2            bdsmatrix_1.3-4            
#> [213] codetools_0.2-16            matrixStats_0.55.0         
#> [215] GO.db_3.10.0                FNN_1.1.3                  
#> [217] gtools_3.8.1                prettyunits_1.0.2          
#> [219] gridBase_0.4-7              GenomeInfoDb_1.22.0        
#> [221] gtable_0.3.0                DBI_1.0.0                  
#> [223] ggpmisc_0.3.5               ggalluvial_0.12.1          
#> [225] highr_0.8                   httr_1.4.1                 
#> [227] KernSmooth_2.23-15          stringi_1.4.3              
#> [229] progress_1.2.2              reshape2_1.4.3             
#> [231] farver_2.0.3                uuid_0.1-4                 
#> [233] viridis_0.5.1               magick_2.2                 
#> [235] coxme_2.2-16                xml2_1.2.2                 
#> [237] rvcheck_0.1.6               boot_1.3-22                
#> [239] lme4_1.1-21                 geneplotter_1.64.0         
#> [241] ggplotify_0.0.4             DESeq2_1.26.0              
#> [243] bit_1.1-14                  ExPosition_2.8.23          
#> [245] ggraph_2.0.0                pkgconfig_2.0.3            
#> [247] mitools_2.4                 tableone_0.12.0

8. Citing MOVICS

MOVICS is a wrapper incorporating algorithms and functions from many other sources. The majority of functions used in MOVICS are incorporated directly from other packages or are drawn from specific published algorithms. Hence, below I provide guidance on which papers should be cited alongside MOVICS when using these functions:

cloud

If you use MOVICS R package in published research, please cite:

  • HAVE NOT BEEN PUBLISHED

Please cite corresponding literature below for multi-omic clustering algorithm if used:

  • CIMLR: Ramazzotti D, Lal A, Wang B, Batzoglou S, Sidow A (2018). Multi-omic tumor data reveal diversity of molecular mechanisms that correlate with survival. Nat Commun, 9(1):4453.
  • iClusterBayes: Mo Q, Shen R, Guo C, Vannucci M, Chan KS, Hilsenbeck SG (2018). A fully Bayesian latent variable model for integrative clustering analysis of multi-type omics data. Biostatistics, 19(1):71-86.
  • Mocluster: Meng C, Helm D, Frejno M, Kuster B (2016). moCluster: Identifying Joint Patterns Across Multiple Omics Data Sets. J Proteome Res, 15(3):755-765.
  • COCA: Hoadley KA, Yau C, Wolf DM, et al (2014). Multiplatform analysis of 12 cancer types reveals molecular classification within and across tissues of origin. Cell, 158(4):929-944.
  • ConsensusClustering: Monti S, Tamayo P, Mesirov J, et al (2003). Consensus Clustering: A Resampling-Based Method for Class Discovery and Visualization of Gene Expression Microarray Data. Mach Learn, 52:91-118.
  • IntNMF: Chalise P, Fridley BL (2017). Integrative clustering of multi-level omic data based on non-negative matrix factorization algorithm. PLoS One, 12(5):e0176278.
  • LRAcluster: Wu D, Wang D, Zhang MQ, Gu J (2015). Fast dimension reduction and integrative clustering of multi-omics data using low-rank approximation: application to cancer molecular classification. BMC Genomics, 16(1):1022.
  • NEMO: Rappoport N, Shamir R (2019). NEMO: cancer subtyping by integration of partial multi-omic data. Bioinformatics, 35(18):3348-3356.
  • PINSPlus: Nguyen H, Shrestha S, Draghici S, Nguyen T (2019). PINSPlus: a tool for tumor subtype discovery in integrated genomic data. Bioinformatics, 35(16):2843-2846.
  • SNF`: Wang B, Mezlini AM, Demir F, et al (2014). Similarity network fusion for aggregating data types on a genomic scale. Nat Methods, 11(3):333-337.

Heatmap generated by MOVICS uses ComplexHeatmap R package, so please cite if any heatmap is created:

  • Gu Z, Eils R, Schlesner M (2016). Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics, 32(18):2847–2849.

If FLAGS is removed when calculating TMB, please cite the following two at the same time, otherwise only the first one:

  • Mayakonda A, Lin D C, Assenov Y, et al. (2018). Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res, 28(11):1747-1756.
  • Shyr C, Tarailo-Graovac M, Gottlieb M, Lee JJ, van Karnebeek C, Wasserman WW. (2014). FLAGS, frequently mutated genes in public exomes. BMC Med Genomics, 7(1): 1-14.

If calculating FGA, please cite:

  • Cerami E, Gao J, Dogrusoz U, et al. (2012). The cBio Cancer Genomics Portal: An Open Platform for Exploring Multidimensional Cancer Genomics Data. Cancer Discov, 2(5):401-404.
  • Gao J, Aksoy B A, Dogrusoz U, et al. (2013). Integrative analysis of complex cancer genomics and clinical profiles using the cBioPortal. Sci Signal, 6(269):pl1-pl1.

Drug sensitivity analysis is based on pRRophetic R package, so please cite:

  • Geeleher P, Cox N, Huang R S (2014). pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PLoS One, 9(9):e107468.
  • Geeleher P, Cox N J, Huang R S (2014). Clinical drug response can be predicted using baseline gene expression levels and in vitro drug sensitivity in cell lines. Genome Biol, 15(3):1-12.

  Please cite corresponding literature below if using differential expression analysis:

  • edgeR:
    • Robinson MD, McCarthy DJ, Smyth GK (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics, 26(1):139-140.
    • McCarthy DJ, Chen Y, Smyth GK (2012). Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation. Nucleic Acids Res. 40(10):4288-4297.
  • DESeq2: Love MI, Huber W, Anders S (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol, 15(12):550-558.
  • limma: Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK (2015). limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res, 43(7):e47.

  Gene set enrichment analysis uses clusterProfiler R package with the following paper:

  • Yu G, Wang L, Han Y, He Q (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS, 16(5):284-287.

  For single sample enrichment analysis, plese cite literatures below appropriately:

  • ssgsea: Barbie, D.A. et al (2009). Systematic RNA interference reveals that oncogenic KRAS-driven cancers require TBK1. Nature, 462(5):108-112.
  • gsva: Hänzelmann, S., Castelo, R. and Guinney, J. (2013). GSVA: Gene set variation analysis for microarray and RNA-Seq data. BMC Bioinformatics, 14(1):7.
  • zscore: Lee, E. et al (2008). Inferring pathway activity toward precise disease classification. PLoS Comp Biol, 4(11):e1000217.
  • plage: Tomfohr, J. et al (2005). Pathway level analysis of gene expression using singular value decomposition. BMC Bioinformatics, 6(1):1-11.

External validation uses nearest template prediction as implanted in CMScaller R package. Please cite:

  • Eide P W, Bruun J, Lothe R A, et al. (2017). CMScaller: an R package for consensus molecular subtyping of colorectal cancer pre-clinical models. Sci Rep, 7(1):1-8.
  • Hoshida, Y. (2010). Nearest Template Prediction: A Single-Sample-Based Flexible Class Prediction with Confidence Assessment. PLoS One, 5(11):e15543.

REFERENCES

  1. Pierre-Jean M, Deleuze J F, Le Floch E, et al (2019). Clustering and variable selection evaluation of 13 unsupervised methods for multi-omics data integration. Brief Bioinformatics.
  2. Rappoport N, Shamir R (2018). Multi-omic and multi-view clustering algorithms: review and cancer benchmark. Nucleic Acids Res, 46(20):10546-10562.
  3. Cancer Genome Atlas Network. Comprehensive molecular portraits of human breast tumours. Nature, 2012, 490(7418):61-78.
  4. Yau C, Esserman L, Moore D H, et al. (2010). A multigene predictor of metastatic outcome in early stage hormone receptor-negative and triple-negative breast cancer. Breast Cancer Res, 5(12):1-15.
  5. Chalise P, Fridley BL (2017). Integrative clustering of multi-level omic data based on non-negative matrix factorization algorithm. PLoS One, 12(5):e0176278.
  6. Tibshirani, R., Walther, G., Hastie, T. (2001). Estimating the number of data clusters via the Gap statistic. J R Stat Soc Series B Stat Methodol, 63(2):411-423.
  7. Le DT, Uram JN, Wang H, et al. (2015). PD-1 Blockade in tumors with mismatch-repair deficiency. N Engl J Med, 372(26):2509–2520.
  8. Davoli T, Uno H, Wooten E C, et al. (2017). Tumor aneuploidy correlates with markers of immune evasion and with reduced response to immunotherapy. Science, 355(6322):261-U75.
  9. Geeleher P, Cox N J, Huang R S (2014). Clinical drug response can be predicted using baseline gene expression levels and in vitro drug sensitivity in cell lines. Genome Biol, 15(3):1-12.
  10. Geeleher P, Cox N, Huang R S (2014). pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PLoS One, 9(9):e107468.
  11. Rand W M. (1971). Objective criteria for the evaluation of clustering methods. J Am Stat Assoc, 66(336):846-850.
  12. Vinh, N. X., Epps, J., Bailey, J. (2009). Information theoretic measures for clusterings comparison. Proceedings of the 26th Annual International Conference on Machine Learning - ICML ’09. p. 1.
  13. Jaccard Distance (Jaccard Index, Jaccard Similarity Coefficient). Dictionary of Bioinformatics and Computational Biology.
  14. Fowlkes E B, Mallows C L. (1983). A method for comparing two hierarchical clusterings. J Am Stat Assoc, 78(383):553-569.
  15. Hoshida, Y. (2010). Nearest Template Prediction: A Single-Sample-Based Flexible Class Prediction with Confidence Assessment. PLoS One, 5(11):e15543.